# 三次样条差值（自然边界条件）

from numpy import zeros, linalg, linspace, diag, ones
from sympy import abc

def one(x, y, a, b, n):
    import matplotlib.pyplot as plt
    h = q = w = e = t = zeros(n)
    u = r = z = zeros(n - 1)
    A = diag(ones(n-1)*2)
    K = g = zeros(n + 1)
    # 计算和，u，r，g
    for i in range(0, n):
        h[i] = x[i + 1] - x[i]
    print("h=", h)
    for i in range(0, n - 1):
        u[i] = h[i] / (h[i] + h[i + 1])
    print("u=", u)
    for i in range(0, n - 1):
        r[i] = 1 - u[1]
    print("r=", r)
    for i in range(1, n):
        g[i] = (6 / (h[i - 1] + h[i])) * ((y[i + 1] - y[i]) / h[i] - (y[i] - y[i - 1]) / h[i - 1])
    g[1] = g[1] - u[0] * a
    g[n - 1] = g[n - 1] - r[n - 2]
    g[0] = a
    g[n] = b
    print("g=", g)
    for i in range(0, n - 1):
        z[i] = g[i + 1]
    print("z=", z)
    # 对线性方程组系数矩阵进行赋值
    for i in range(0, n - 2):
        A[i + 1][i] = u[i + 1]
    for i in range(1, n - 1):
        A[i - 1][i] = r[i - 1]
    print("A=", A)
    # 求解线性方程组
    M = linalg.solve(A, z)
    for i in range(1, n):
        K[i] = M[i - 1]
    K[0] = a
    K[n] = b
    print("K=", K)
    for i in range(0, n):
        q[i] = (K[i + 1] - K[i]) / 6 * h[i]
        w[i] = (K[i] * x[i + 1] - K[i + 1] * x[i]) / 2 * h[i]
        e[i] = (3 * K[i + 1] * pow(x[i], 2) - 3 * K[i] * pow(x[i + 1], 2) - 6 * y[i] + K[i]
                * pow(h[i], 2) + 6 * y[i + 1] - K[i + 1] * pow(h[i], 2)) / 6 * h[i]
        t[i] = (K[i] * pow(x[i + 1], 3) - K[i + 1] * pow(x[i], 3) + 6 * y[i] * x[i + 1] - K[i]
                * pow(h[i], 2) * x[i + 1] - 6 * y[i + 1] * x[i] + K[i + 1] * pow(h[i], 2) * x[i]) / 6 * h[i]
    print("q=", q)
    print("w=", w)
    print("e=", e)
    print("t=", t)

    # 对区间进行划分
    x1 = linspace(0, 0.4, 100)
    x2 = linspace(0.4, 0.8, 100)
    x3 = linspace(0.8, 1.2, 100)

    # 分段函数定义
    y1 = (q[0] * pow(x1, 3) + w[0] * pow(x1, 2) + e[0] * x1 + t[0])
    y2 = (q[1] * pow(x2, 3) + w[1] * pow(x2, 2) + e[1] * x2 + t[1])
    y3 = (q[2] * pow(x3, 3) + w[2] * pow(x3, 2) + e[2] * x3 + t[2])
    print("S(x1)=" + str(q[0]) + "x^3+" + str(w[0]) + "x^2+" + str(e[0]), "x+" + str(t[0]))
    print("S(x2)=" + str(q[1]) + "x^3+" + str(w[1]) + "x^2+" + str(e[1]), "x+" + str(t[1]))
    print("S(x3)=" + str(q[2]) + "x^3+" + str(w[2]) + "x^2+" + str(e[2]), "x+" + str(t[0]))

    # 画出图像
    plt.plot(x1, y1)
    plt.plot(x2, y2)
    plt.plot(x3, y3)
    plt.show()


# 主函数部分
x = linspace(0.2, 1.0, 5)
y = [0.98, 0.92, 0.81, 0.64, 0.38]
n = len(x) - 1
one(x, y, 1, -1, n)

datax=[0.2+0.08*i for i in [0,1,11,10]]

datay=[]
for k in datax:
    if 0<=k<0.4 :
        y11=
